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We present and test a new method to compute the hadronic vacuum polarization function in 
lattice simulations. This can then be used, e.g., to determine the leading hadronic contribution to 
the anomalous magnetic moment of the muon. The method is based on computing susceptibilities 
with respect to external electromagnetic plane wave fields and allows for a precision determination of 
both the connected and the disconnected contributions to the vacuum polarization. We demonstrate 
that the statistical errors obtained with our method are much smaller than those quoted in previous 
lattice studies, primarily due to a very effective suppression of the errors of the disconnected terms. 

These turn out to vanish within small errors, enabling us to quote an upper limit. We also comment 
on the accuracy of the vacuum polarization function determined from present experimental i?-ratio 
data. 
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I. INTRODUCTION 

The most precise measurement of the anomalous 
magnetic moment of the muon, obtained by E821 at 
Brookhaven [1], differs by more than three standard de¬ 
viations from the theoretical expectation. At present, 
the uncertainties on the theory and on the experimental 
sides are of similar sizes. For recent reviews and analy¬ 
ses, see, e.g., refs. [2-6]. With the planned E989 exper¬ 
iment at Fermilab [2] and E34 at J-PARC [7], it is of 
utmost importance to increase the precision of the stan¬ 
dard model prediction in line with the expected experi¬ 
mental improvement by a factor of about five [1, 2, 7]. If 
the discrepancy persisted at this even higher level of accu¬ 
racy, this should help to pin down any particular beyond- 
the-standard-model scenario and constrain the parame¬ 
ters of new interactions. With an impressive QED five- 
loop evaluation [8] available, the theoretical uncertainty 
is dominated by non-perturbative effects and, in partic¬ 
ular, by the leading hadronic contribution to the elec¬ 
tromagnetic vacuum polarization tensor, with the second 
biggest source of uncertainty being the hadronic light-by¬ 
light scattering contribution. The hadronic contribution 
to the vacuum polarization tensor is also important in 
view of the running of the electromagnetic fine structure 
constant and of the Weinberg weak mixing angle [4, 9-12] 
from low to high scales. 

The standard method [13, 14] employed in lattice cal¬ 
culations of the leading hadronic contribution to the 
anomalous magnetic moment — 2)/2 of a charged 

lepton £ S {e,p-, r} consists of computing the renor¬ 
malized vacuum polarization function and inserting this 
into the leading-order QED formula. The hadronic vac¬ 
uum polarization tensor, which is the main object of this 
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study, is defined as 

= J dt: (>(a;)>(0)) (1) 

= {PuPy - ^urP^) n(p^), 

where 

Qu - , J j 9s _ /r,\ 

Ju = —u'y^u -\ -H-H- (2) 

e e 

denotes the quark electromagnetic current in position 
space and qu/e = 2/3, 9 d/e = 9 s/e = —1/3 are the 
fractional quark charges. Due to electromagnetic current 
conservation, 11,^^ is transverse and can be parameterized 
in terms of a single vacuum polarization function 11 (p^), 
where we employ Euclidean spacetime conventions, i.e. 
the spacelike > 0 correspond to virtualities, n(p^) 
undergoes additive renormalization but the renormalized 
combination 

BrCp") = n(p2) - n(0) (3) 

is ultraviolet finite. 

It turns out that the leading hadronic contribution 
^had,LO .j-Q |.]^g anomalous magnetic moment of the muon 
(see the definition equation (4) below) depends most 
strongly on nR(p^) at relatively small argument values. 
Since small momenta correspond to large Euclidean dis¬ 
tances, naively implementing equation (1) results in a 
bad signal over noise ratio in this region. This becomes 
even worse for calculations of the quark-line disconnected 
contributions, which therefore have been neglected in al¬ 
most all previous lattice studies. Where these were taken 
into account [15-17], they dominated the statistical er¬ 
ror. Another problem of many past lattice attempts is a 
conceptual one: 11(0) often is extrapolated from n(p^) at 
p^ > 0 and the parametrization used constitutes a source 
of systematic uncertainty that is difficult to estimate. 

Here we propose methods that address both of the 
above issues. The vacuum polarization at p^ = 0 is shown 
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to be equal to the bare magnetic susceptibility of the 
system, which can be determined independently on the 
lattice. We investigate different methods to achieve this, 
giving consistent results. We also discuss how this quan¬ 
tity diverges as a function of the lattice spacing towards 
the continuum limit. 

Most importantly, we introduce a new method for com¬ 
puting both the connected and the disconnected con¬ 
tributions to the hadronic vacuum polarization function 
with unprecedented precision, in particular at small mo¬ 
menta. This consists of calculating n(p^) at > 0 
through the response of the system to oscillatory back¬ 
ground electromagnetic fields. The new method is sim¬ 
ilar in spirit to employing momentum sources [18, 19], 
allowing us to spend more effort on the low-p^ points, 
thereby increasing their precision, without wasting re¬ 
sources on large momenta where nR(p^) can easily be 
obtained within small relative errors, with a much smaller 
impact on the predicted value of 

The methods are tested on Nf = 2 + 1 staggered en¬ 
sembles at the physical point, neglecting QED effects 
on the quark propagation which are of a higher order 
in the fine-structure constant a. In this situation, due 
to J2f^{uds}^f ~ disconnected contributions van¬ 
ish for nis = rUud but need to be taken into account 
for rus > nriud, which we do. Since we neglect charm 
quark effects, we have to restrict ourselves to < m^. 
At high momenta our results can, however, be combined 
with measurements of the i?-ratio as well as with per¬ 
turbation theory: the non-singlet and singlet QCD con¬ 
tributions to the Adler function have been calculated in 
massless QCD to 0(af) in the strong coupling constant 
in refs. [20] and [21], respectively. 

This article is organized as follows. In section II we 
review previous calculational strategies, followed by sec¬ 
tion III, where we introduce our background field method 
and link this to magnetic susceptibilities. We also dis¬ 
cuss renormalization issues and comment on relations 
between the Adler function and the entropy density at 
high temperatures. Finally, in section IV we present the 
simulation setup and first results, before we conclude. 
The equivalence between magnetic susceptibilities and 
the vacuum polarization is demonstrated in appendix A, 
and the details of our numerical implementation are dis¬ 
cussed in appendices B and C. 


II. SUMMARY OF PREVIOUSLY EMPLOYED 
METHODS 


The leading hadronic contribution to the anomalous 
magnetic moment is given as [13, 22] 

^had.LO ^ ^^2 / dp2A:£(p^)nR(p2) , (4) 

^0 

where the perturbative kernel function is defined as 


Keip^) 


mjp^Zi{p^)^ [l - p^Zi,{p^)] 
I -I- m?^p^Z{p^Y 


(5) 


with 


Zi{p^) 


^/l + - I 

2m^ 


( 6 ) 


The renormalized hadronic vacuum polarization function 
is defined in equations (1) and (3) above. Note that the 
above expressions are valid to leading-order in terms of 
the QED fine-structure constant a = e^/(47r) « 1/137, 
i.e. to which, at this order, can be pulled out of 

the integral. 

In the limit of small momenta, where nn(p^) cx p^, 
the argument of the integral has its maximum at pg « 
(•\/5 — 2)rrig. For the muon with « 0.105 GeV this 
implies pg « 0.0026 GeV^: an enormous volume would 
be necessary to resolve this momentum region, at least 
without the use of twisted boundary conditions [23, 24], 
since tt/pq « 27r/m^ « 12 fm. Fortunately, the inte¬ 
gral as a whole turns out to be dominated by somewhat 
higher momenta: it still picks up about 50% of its value 
from momenta larger than 10 pg. The predicted value 
of strongly depends on nR(p^) at these still rel¬ 

atively small momenta p^ ~ 0.03 GeV^. This is nicely 
illustrated, e.g., in ref. [25], in figure 3 of ref. [24] and in 
figure 1 of ref. [17]. 


A. Information from experiment 


The hadronic polarization tensor (and also the lead¬ 
ing hadronic contribution to the lepton anomalous mag¬ 
netic moments [22]) can be obtained by analytic contin¬ 
uation of the i?-ratio of the total cross section a{e'^e~ —)■ 
hadrons) over the tree-level e“''e“ —>■ expectation 

(see, e.g., refs. [26, 27]): 


nR(p") 


127r2 


S(s+p2)- 


(7) 


i?-ratio measurements [4, 5] can in principle be aug¬ 
mented by other experimental data, including r-decays 
into final states containing tt+tt”, see, e.g., ref. [28]. 

In figure 1 we show the so-determined renormalized 
vacuum polarization as a function of p^ [29] . The present 
relative precision of Hr is 0.64% at p^ = 0.025 GeV^, in¬ 
creasing to 0.74% at p2 = 0.6GeV^ [5, 29]. Achieving 
a statistical error below 1% around p^ = 0.03 GeV^ al¬ 
ready constitutes an enormous challenge for present-day 
lattice determinations, and such results still need to be 
extrapolated to the infinite volume and continuum limits 
and, often, to physical quark masses. In principle, lat¬ 
tice data at large p^ values - where discretization errors 
are enhanced - can be substituted by results from the R- 
ratio. Such a combined strategy may prove optimal for 
an accurate determination of once sufficiently 

precise lattice results become available. 
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FIG. 1. The renormalized vacuum polarization deter¬ 
mined from the experimental iZ-ratio [29], performing the in¬ 
tegral (7) up to s = Smax = (2GeV)^, where three quark 
flavors are active. Also indicated is the result of the in¬ 
tegral supplemented by three-flavor perturbation theory for 
s > (2GeV)^ 

B. Lattice determinations of 11(0) 

In the past, two strategies have been used to obtain 
the zero-momentum subtraction 11(0). One possibility 
are hts of n(p^) data, e.g., to pole parametrizations, as¬ 
suming vector dominance [14, 24, 30-32], which is also 
suggested to be the dominant contribution by chiral per¬ 
turbation theory [30]. Extending the ht region towards 
large momenta, such pole ansatze have also been com¬ 
bined with polynomial parametrizations [15, 17, 30, 33], 
motivated by perturbation theory. Another popular and 
less model-dependent way to obtain the normalization is 
through Fade approximants [24, 25, 34, 35]. 

As an alternative, one can compute derivatives of 
n^iy(p) from its definition in terms of the continuum 
Fourier transformation (1). Then the divergent contri¬ 
bution that needs to be subtracted from n(p^) can, e.g., 
be obtained via 

n(0) = ~ 9 )^7wn;yi/(p) ifJ- 7^ i') 

p=0 

= \ J (>(x)>(0)) = i Jdtf Git ), (8) 

where no summation over v is implied and in the last 
step we identified /r with the time-direction, to emphasize 
the correspondence to the second t-moment of a zero- 
momentum projected two-point function 

Git) = Jd\ iji (r, t)j, (0)). (9) 

This method was used, e.g., in refs. [36, 37], to obtain 
this subtraction. 

Finally, in ref. [38] the expansion of the two-point 
current-current correlation function in powers of is 


carried out already on the level of quark propaga¬ 
tors. This enables the direct computation of 11(0) = 
5^ni2/(i9pi9p2)|p=0) without relying on a continuum for¬ 
mula. However, this comes at the price of computing 
the expectation value of an operator involving up to four 
fermion matrix inversions, without even considering dis¬ 
connected contributions. 

C. Lattice determinations of n(p^), IIr(p^) or 
moments thereof 

The lattice vector Ward-Takahashi identity reads 
= 0 and therefore [13, 14, 39] 

^pyip^) = {Pp.Pu - n(p^), (10) 

where Pp = (2/a) sin(ap^/2). This change that affects 
n(p^) at high momenta has been implemented in almost 
all lattice studies, as well as a modified phase i—7 

gip(a:+a/i/ 2 -a!y/ 2 ) Fourier sum for p ^ V. 

Most lattice evaluations use what we will refer to be¬ 
low as the conventional method. This amounts to di¬ 
rectly computing the lattice version of equation (1), see, 
e.g., refs. [13-15, 30-33, 37]. In some cases, lower mo¬ 
menta were made accessible by the use of twisted bound¬ 
ary conditions [24, 40, 41]. Very recently, another in¬ 
teresting method, stochastically averaging over different 
twists, has been suggested [42] that reduces finite vol¬ 
ume effects and allows to realize very small momenta. 
The main problem of modifying the fermionic boundary 
conditions is that this cannot easily be extended to in¬ 
corporate quark-line disconnected contributions. 

Obviously, equation (4) can be Taylor expanded in 
powers of p^ and the coefficients can be related to those 
of the corresponding expansion of nR(p^). Generalizing 
equation (8) above, the first and higher order derivatives 
of with respect to p^ can be obtained, computing 
t^-moments of two-point zero-momentum (spatial) pro¬ 
jected current-current correlators. This was explored 
within ref. [43] and carried out for the first few moments 
of the connected strange and charm quark contributions 
to in ref. [44]. In ref. [27] the anomalous magnetic 

moment was directly related to the zero momentum pro¬ 
jected current-current two-point function. This approach 
was then employed, e.g., in refs. [16, 36, 41]. 

So far, disconnected contributions have been included 
in very few lattice studies [15-17]. While their effect 
seems to be small, the associated statistical error ex¬ 
ceeds that of the connected terms. Here we will find 
that this need not be the case. There exist theoreti¬ 
cal expectations regarding the size of flavor singlet con¬ 
tributions: exploiting the fact that m^,m^ > rrip, it 
was demonstrated [16] that the ratio of the disconnected 
contribution over the total momentum-projected current- 
current two-point function Git), defined in equation (9), 
approaches the value —1/9, in the limit of large Euclidean 
times for Nf = 2-1-1 quark flavors. This ratio will, 
however, not automatically propagate into nR(p^) that 
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depends on G{t) at all times t. Next-to-leading order 
chiral perturbation theory arguments show the discon¬ 
nected contribution to also account for —1/9 of the to¬ 
tal nH(p^) [45]. However, this observation builds on the 
fact that the correlator of the iso-singlet vector current 
ujfiU -I- is momentum-independent to this order of 
chiral perturbation theory — which we found is not at all 
satisfied by the lattice data. Thus, direct computation of 
the disconnected terms cannot be avoided in a systematic 
study. Our numerical results will shed light onto the size 
of the disconnected contribution at low . 

III. VACUUM POLARIZATION FROM 
SUSCEPTIBILITIES 

A. The method 


since in Euclidean spacetime at zero temperature 
the indices can be relabelled at will. 

To find the background held corresponding to n^j,(p), 
we dehne the magnetic helds 

HP (x) = B siii{px) Bs , B° = ile 3 , (13) 

pointing in the third spatial direction. While is an os¬ 
cillatory magnetic held with oscillation frequency p, B*^ 
is a homogeneous background. The corresponding sus¬ 
ceptibilities are obtained as the second derivatives of the 
free energy density with respect to the amplitude of the 
magnetic held: 


d^mp] 

d{eBY 


B=0 


(14) 


The photon vacuum polarization tensor (1) can also be 
interpreted as a momentum space current-current corre¬ 
lation function 

^ , (11) 

where V 4 denotes the four-dimensional volume of the sys¬ 
tem and jfi is the Fourier transform of the electromag¬ 
netic current dehned in equation ( 2 ): 

J/^b) = j d'^xe^P^j^ix). ( 12 ) 

Depending on the lattice dehnition of j/j, the polarization 
tensor ( 11 ) may or may not renormalize multiplicatively 
with Zy. Here, we work with a conserved current, i.e. 
Zy = 1. 

In the following we will relate the vacuum polariza¬ 
tion to the leading response of the free energy density / 
of the system to background electromagnetic helds. To 
illustrate the relation between the two objects on a qual¬ 
itative level, it is instructive to represent the vacuum 
polarization tensor by the diagram 

p V 

VWVA 


where a momentum p hows in and out of the photon legs. 
Here, the gray blob indicates all possible closed loops 
formed by quark and gluon propagators — i.e. the per¬ 
turbative expression for the free energy density /. The 
legs may be thought of as photons corresponding to a 
background electromagnetic held with momentum p. 
Pulling out these legs is achieved by taking appropri¬ 
ate derivatives of / with respect to the background held. 
While background electric helds turn the Euclidean QCD 
action complex and are thus problematic in lattice simu¬ 
lations, background magnetic helds can be realized with¬ 
out complications. Employing the latter gives access to 
the spatial components H^ and hence to all components 


These susceptibilities are normalized by the square of the 
elementary charge e > 0 to ensure that only the renor¬ 
malization group-invariant combination eB appears in 
the dehnitions. Note that Xp can be evaluated on gauge 
ensembles generated at H = 0 . 

The explicit calculation in appendix A shows that 

2xp = n(p2), = n(0). (15) 

These relations form a new representation of the vacuum 
polarization function in terms of susceptibilities with re¬ 
spect to the magnetic helds dehned in equation (13) and 
are the main result of this article. 

Unlike the conventional method, where the polariza¬ 
tion function is extracted from the same set of posi¬ 
tion space current-current correlators for all momenta, 
equation (15) gives access to n(p^) at one single lattice 
momentum p. While this certainly increases the costs 
of calculating H over a large range of momenta, it also 
allows for a better signal-to-noise ratio within momen¬ 
tum regions of particular interest. As argued above, 
for the determination of the hadronic contribution to 
the muon anomalous magnetic moment low mo¬ 

menta p^ ^ 0.03 GeV^ are much more important than 
the high-p region. While {j^{x)j^(d)) mixes informa¬ 
tion about all allowed values of p, here such a mixing 
is avoided. 

Just as the vacuum polarization tensor, Xp and xo can 
also be separated into connected and disconnected con¬ 
tributions. We will demonstrate in section IV below that, 
using this new approach, an unprecedented accuracy can 
be achieved for both the connected and the disconnected 
contributions to the vacuum polarization function, al¬ 
ready at moderate computational costs. An additional 
advantage of the method is that it gives direct access to 

n(o). 

To summarize, to arrive at a prediction for it is 

desirable to improve the accuracy in the low-p region and 
to calculate n(0) independently. The method we propose 
accomplishes both of these requirements. 
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B. Renormalization 

Before presenting the details of the implementation 
and our numerical results, it is instructive to discuss the 
renormalization properties of xo in more detail. Equa¬ 
tion (15) reveals that the homogeneous susceptibility is 
additively divergent, just as n(0). To see where this di¬ 
vergence comes from, let us consider the multiplicative 
renormalization of the background magnetic field (and 
the corresponding renormalization of the electric charge), 

e^ = Z~^el, B'^ = ZeBl, eB = CrBr, (16) 

with the renormalization factor 

Zf, = l + 2biel log(/ra), (17) 

where a is the lattice spacing (inverse of the regulator) 
and fi the renormalization scale. Notice that since the 
magnetic field is external and has no dynamics, only the 
lowest-order QED /3-function coefficient — denoted as bi 
— appears in Z,, [46-48]. 

The total free energy density of the system is the 
sum of / and the energy i?^/2 of the magnetic field. Since 
varying the background field should not change the ul¬ 
traviolet properties of the system, must be free of 
B-dependent divergences. This implies that the diver¬ 
gence of the pure magnetic energy 

^ ^ -b 6i {eBf \og{fia) (18) 

is exactly cancelled by an analogous divergence of /. 
Plugging this divergence into the definition (14), we ob¬ 
tain 


Xo = 26i(a) log(/xa). (19) 


The renormalization scale /i is fixed by the requirement 
that there should be no finite quadratic terms in 
other than i?^/2 [46]. Let us emphasize that bi is the 
lowest-order coefficient of the QED /3-function, however, 
with all QCD corrections taken into account. To high¬ 
light this, we explicitly indicate the dependence of bi on 
the lattice spacing. Perturbatively, this reads [49] 


bi{a) 


f — U,d,S 


1-b 


aHa) 

47r^ 



( 20 ) 


where g^{a) is the QCD coupling. Equation (19) allows 
to connect lattice results for xo to perturbation theory, 
once the lattice spacing is small enough, cf. ref. [48]. 


is defined as the logarithmic derivative of the polarization 
function with respect to the squared momentum [26]: 


D{p^) = 127r2 


an(p2) 
aiog(p2) ■ 


( 21 ) 


Let us consider QCD at a high temperature T, which 
exceeds all other dimensionful scales in the system. In 
this limit, the argument of 11 is set by a thermal scale 
^th = 27rT, leading to the correspondence n(/rjjj) o 
Xo(T^)- (The susceptibility at high temperatures indeed 
only depends on [48].) For the Adler function, this 
implies the relation 




127r" 


dxo 

aiogr2 


= 


d'^i 


d{eBy 


( 22 ) 


B=0 


where in the second step we used the definition of the 
entropy density s = —df/dT. Equation (22) reveals 
that the leading dependence of the entropy density on 
the magnetic field at high temperatures is fixed by the 
Adler function, i.e. by perturbative QED physics. Re¬ 
peating the above argument with T replaced by a chem¬ 
ical potential /i (or by an isospin chemical potential /ii) 
gives an analogous relation for the quark number den¬ 
sity n = —df/dfi at high p, (or for the isospin density 
ni = —df/dfii at high ^i). We believe these are highly 
non-trivial findings. 


IV. SIMULATION DETAILS AND NUMERICAL 
RESULTS 

We employ the Nf = 2 -|- 1 staggered lattice ensem¬ 
bles [50, 51] generated at physical pion and kaon masses. 
Each ensemble — summarized in table I — consists of a 
hundred to a few hundred effectively statistically decorre- 
lated configurations. Details of the simulation algorithm 
and of the lattice setup can be found in refs. [50, 52, 53]. 


TABLE I. Lattice ensembles investigated; the largest lattice 
spacing reads ao = 0.29 fm. 


w 

Nt 

P 

a [fm] 

log(a/ao) 

24 

32 

3.45 

0.290 

0 

24 

32 

3.55 

0.216 

-0.295 

32 

48 

3.67 

0.153 

-0.636 

40 

48 

3.75 

0.125 

-0.843 

40 

48 

3.85 

0.099 

-1.078 


C. Implication for hot or dense QCD 

As a side-remark, we mention that the correspon¬ 
dence (15) can be generalized to high temperatures. In 
this case it results in a relation between the entropy den¬ 
sity and the perturbative Adler function [48]. The latter 


A. Oscillatory susceptibilities 

First we discuss results on the susceptibilities Xp = 
n(p^)/2 with respect to the oscillatory backgrounds. 
These are determined via the noisy estimator tech¬ 
nique described in appendix B. A typical set of low- 
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FIG. 2. The low-momentum region of the oscillatory sus¬ 
ceptibilities as measured on the 24® x 32 configurations at 
P = 3.45. The curves correspond to polynomial- and Pade- 
type extrapolations of 2xp to p = 0. The direct determination 
Xo is shifted horizontally to the left for better visibility. Also 
included are results obtained using random wall sources, dis¬ 
placed horizontally to the right. 


momentum results is shown in figure 2. The data in¬ 
clude both the connected and the disconnected contri¬ 
butions to n(p^). The figure also includes results ob¬ 
tained via the conventional method, however, employing 
stochastic wall sources (for our numerical implementa¬ 
tion, see appendix C). The comparison reveals full agree¬ 
ment between the two approaches. The statistical er¬ 
ror of the random wall data increases towards small mo¬ 
menta, whereas it remains tiny even for the lowest non¬ 
vanishing p^-value shown for the oscillatory susceptibil¬ 
ities. Note that the number of inversions employed to 
obtain the data point at the lowest momentum was the 
same, iVinv = 3000, for both approaches. 

In most previous lattice studies, 11(0) was obtained 
by extrapolating Tl{p^) to zero. Some possible extrap¬ 
olations, employing polynomials or Fade approximants, 
fitted over various ranges in are included in the fig¬ 
ure. These fits are also compared to the direct determi¬ 
nations via the homogeneous susceptibility xo (see sec¬ 
tion IV B below) and via the zero-momentum projected 
current-current correlation function G(t) according to 
equation (8), again obtained using random wall sources. 
Within their scatter, at = 0 the extrapolations agree 
with the direct determinations. We remark that increas¬ 
ing the precision for the lowest few momenta stabilizes 
such extrapolations tremendously. 


B. Homogeneous susceptibility and renormalized 
vacuum polarization 

The susceptibility xo with respect to a homogeneous 
background is of interest for QCD thermodynamics in 
magnetic fields and has been the subject of detailed stud¬ 


ies in the past few years. The determination of xo is 
considerably more complicated than that of Xp due to 
the quantization of the magnetic flux $. On the one 
hand, oscillatory magnetic fields have zero flux and can 
be varied continuously, allowing for a direct differentia¬ 
tion with respect to B. On the other hand, homogeneous 
fields have nonzero flux. Therefore, such a differentiation 
cannot be carried out to determine xo: see appendix B. 
Several approaches, summarized in refs. [48, 54], have 
been developed recently to overcome this problem. Here 
we compare results obtained using the finite difference 
method [55], the generalized integral method [48] and 
the half-half method [56]. The former two approaches 
are based on simulations at non-zero magnetic flux val¬ 
ues, numerically differentiating the results with respect 
to <i>. The half-half method involves calculating expecta¬ 
tion values directly at i? = 0, employing a setup where 
the magnetic field is positive in one half and negative in 
the other half of the lattice. In this case, since the total 
flux is zero, a direct differentiation with respect to the 
amplitude is possible. However, the discontinuity of the 
magnetic field turns out to dramatically enhance finite 
volume effects in xo; see below. ^ 

In figure 3, we compare all three approaches. The re¬ 
sults from the generalized integral method and from the 
finite difference approach are taken from refs. [58, 59] 
while the half-half results are new. Not all lattice spac- 
ings are covered by all the methods. While the results of 


-0.06 

-0.07 

-0.08 

-0.09 


FIG. 3. Magnetic susceptibility with respect to a homoge¬ 
neous background as a function of the logarithm of the lat¬ 
tice spacing (ao = 0.29 fm), using three different approaches 
(the generalized integral method [48], the finite difference 
method [58, 59] and data generated in this study using the 
half-half method [56]). Also included are a comparison to 
0{g^) perturbation theory and a parametrization via a ratio¬ 
nal ansatz. 


Zi I I I I I I I I I I I I I I I I I I I I I I r 


I gen. int. method [1406.0269] 
$ finite diff. [1310.8656] e 
^ half-half method . /■' 


- IE 


I- 1,-r I 


."5 


- - ©(g^) scaling 
■ O(g^) X rational 


_L 


_L 


_L 


_L 


■1 -0.8 -0.6 -0.4 -0.2 0 

log(o/o„) 


^ These finite volume effects cancel to a large extent in the differ¬ 
ence xo(^) ~ XoC^ = 0) [57], which is relevant for QCD thermo¬ 
dynamics in background magnetic fields. 
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the generalized integral method^ and of the finite differ¬ 
ence approach are consistent with each other, the half¬ 
half approach consistently underestimates the magnitude 
of the susceptibility. The difference between that ap¬ 
proach on the one hand and the other two methods on 
the other hand is found to be as large as 10% and re¬ 
duces only very slowly with increasing lattice volumes.^ 
Altogether, we conclude that the half-half method is in¬ 
sufficient for our purposes and discard it in the following. 

Perturbation theory predicts the dependence of xo on 
the lattice spacing, see equations (19) and (20). In fig¬ 
ure 3 the data are plotted against log(a/ao) to verify the 
expected logarithmic divergence. We include the lead¬ 
ing 0{g'^) QCD correction to the lowest-order QED (3- 
function coefficient bi. The renormalization scale g is 
fitted to match the lattice results (dashed green line). In 
addition, we multiply the resulting curve by a rational 
function that approaches unity as a —> 0 (solid yellow er¬ 
ror band). This band defines the homogeneous magnetic 
susceptibility Xo(a)) as shown for one lattice spacing in 
the very left of figure 2. The resulting renormalization 
scale reads g = 0.123(8) GeV, consistent with our deter¬ 
mination in ref. [48] . 

The n(p^) results are shown for all five ensembles of 
table I in figure 4, where 11(0) = xo with the suscepti¬ 
bility Xo determined as detailed above. Notice that the 
statistical uncertainties (again, both connected and dis¬ 
connected terms are taken into account) within our win¬ 
dow of lattice spacings remain at the sub-percent level 
for > 0 and are about one percent for p = 0. Taking 
into account the statistical errors of n(p^) and of the in¬ 
dependently determined 11(0), the renormalized vacuum 
polarization (3) is plotted in figure 5 for the whole mo¬ 
mentum region under consideration. For orientation we 
also show the three flavor perturbation theory result for 

> 2 GeV^, where we truncate the formulae of refs. [20] 
and [21] at 0{al). The perturbative curve is only defined 
up to an overall constant shift, which we adjust by match¬ 
ing to a continuum extrapolation around p^ = 2GeV^. 
It is clear from the figure that — as one would expect — 
lattice spacing effects become more prominent towards 
high momenta. In addition, the vacuum polarization ob¬ 
tained from the experimental i?-ratio (cf. the blue points 
in figure 1) is included in figure 5. 

Having obtained the renormalized hadronic vacuum 
polarization, we can use equations (4) to (6) [13, 22] to 
predict its contribution to the muon anomalous magnetic 
moment. Choosing a third-order spline interpolation, we 


^ Here we compare data obtained on Nt > Ns zero-temperature 
lattices. On the configurations of ref. [48] at finite (but low) 
temperatures, xo was found to have slightly smaller absolute 
values for fine lattices of table I (/3 > 3.67). 

® The comparison between the half-half method and the general¬ 
ized integral method on our coarsest lattice, already presented 
in ref. [48], has been updated by increasing the statistics and the 
number of noisy estimators to reveal the significant difference 
visible in figure 3. 



p2 (GeV2) 

FIG. 4. Vacuum polarization via magnetic susceptibilities in 
the low-momentnm region. The data include both connected 
and disconnected contributions. 



FIG. 5. Subtracted vacuum polarization with independent 
determinations of n(p^) and 11(0). The data include both 
connected and disconnected contributions. The solid red line 
indicates the experimental result (cf. figure 1) and the dotted 
line the three-loop perturbative prediction (see the text). 


obtain values in the range = (4... 5) • 10“® and 

an upward trend towards the continuum limit. This is 
encouraging as the i?-ratio predictions of refs. [5] and [4] 
for the Nf = 2 -|- 1 -I- 1 flavor theory read = 

6.923(42) • 10-® and = 6.949(43) • 10"®, respec¬ 

tively. However, given that the present lattices are rather 
coarse (0.1 fm < a < 0.3fm), we do not yet attempt a 
full-fledged continuum limit extrapolation. (Note that at 
these lattice spacings, the taste splitting of the staggered 
pion multiplet is still sizeable [53]. Thus, large lattice 
artefacts originating from the heavier pion states are not 
unexpected, since is highly sensitive to the pseu¬ 

doscalar masses.) 








C. Statistical accuracy and disconnected 
contributions 

Next, we perform a quantitative comparison between 
the oscillatory susceptibility method, the conventional 
approach with random wall sources and that with point 
sources. We demonstrate that the statistical error of 
n(p^) can be pushed well below that of existing stud¬ 
ies in the literature - even with the disconnected terms 
taken into account. 

We calculated n(p^) using all three methods on 120 
configurations from the (3 = 3.45 ensemble for a single 
momentum = 0.03 GeV^ using an increased number 
of sources. Figure 6 shows the statistical error as a func¬ 
tion of the number of inversions Wnv The details of our 
implementation can be found in appendices B and C. As 
visible in the figure, the oscillatory susceptibility method 
allows to save 50 — 60% of the computational effort with 
respect to the random wall approach. This difference 
mainly comes from the disconnected contributions, which 
can be calculated very accurately via susceptibilities. In 
fact, the statistical error in this approach is dominated 
by the connected contribution,"^ as is also visible in the 
figure. As expected, the conventional method with point 
sources is not applicable for the determination of the dis¬ 
connected terms. Obviously, it is favorable in terms of 
the total computer time spent to increase the number 
of configurations instead of the number of inversions per 
configuration. We remark that the total number of exact 
inversions necessary to achieve a given error can be con¬ 
siderably reduced by methods like the hopping parameter 
expansion [60, 61], truncated eigenmode substitution [62- 
64], the truncated solver method [65-67] and, in the case 
of Wilson-like fermions, employing spin-explicit stochas¬ 
tic sources [68-70]. 

Finally, we discuss the disconnected contribution 
in more detail. A particular feature of is that it re¬ 
quires no additive renormalization. To see this, note that 
n‘^'s(0) vanishes in the perturbative continuum limit, 
since it is of order g^{a) in the strong coupling [21], 
which dampens the logarithmic divergence and results 
in n‘^‘®(0) to fall off as 1/ log^(o) for a —)■ 0. In our three- 
flavor case the disconnected term even vanishes identi¬ 
cally in perturbation theory due to J2f=u d s^f ~ once 
quark masses can be neglected, i.e. a~^ m^. Based 

on this observation, in figure 7 we plot the unsubtracted 
disconnected vacuum polarization for all our lattice spac- 
ings. (The number of inversions was Wnv = 800 for each 
momentum, with the exception of the left-most point.) 


^ To see why this is the case, note that the number of estimates 
increases quadratically with Afjnv for the disconnected terms but 
only linearly for the connected ones, see the discussion in ap¬ 
pendix B. Therefore, the error on the latter eventually overtakes 
that of the former, before both show the expected asymptotic 
(T^ ~ ci(l -f C 2 /Afinv) fall-off. The inherent gauge noise ci can 
only be reduced by increasing the number of configurations. 



FIG. 6. Statistical error of the total (connected plus dis¬ 
connected) n(p^ = 0.03GeV^) as a function of the number 
of inversions. Compared are the results obtained from oscil¬ 
latory susceptibilities, using point sources and random wall 
sources. In addition, the error of the connected oscillatory 
susceptibility alone is shown. Note the logarithmic scale. 


Overall, is consistent with zero, where the two points 
that deviate by more than two standard deviations from 
this assumption are statistically expected and no system¬ 
atic dependence on the lattice spacing or on the volume 
is apparent. With the exception of three outliers with 
large error bars, all central values are below 2 • 10“^ in 
magnitude. 


0.002 


0.001 

■■6 

c 

0 


FIG. 7. Disconnected contribution to n(p^) as a function of 
for our five lattice spacings. 

Using all available estimators (M^v = 20 000) for the 
P = 3.45 ensemble at p^ = 0.03 GeV^, our most accurate 
determinations for the unsubtracted and the subtracted 
vacuum polarizations read 

p^ = 0.03 GeV^ : U = -0.058362(117), 

= -b0.00002I(026), (23) 

Ur = -1-0.002355(198). 
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Here, n(p^) and were measured using the oscil¬ 

latory susceptibility method. (We highlight again that 
the error of H'^*® is much smaller than that of the total 
n.) The vacuum polarization at zero momentum was ob¬ 
tained via random wall sources. Based on the discussion 
above about the vanishing of n'^“(0) in the continuum 
limit, only the connected part of n(0) is necessary for 
the subtraction. The relative error of the so-obtained 
Hr at this momentum is 8%, and is dominated by the 
error of n(0). Clearly, towards higher p^, where the mag¬ 
nitude of n(p^) increases, the relative error on Hj^ rapidly 
decreases. 


V. SUMMARY 

We developed a new approach to determine the 
hadronic vacuum polarization n(p^) on the lattice. It 
is based on calculating magnetic susceptibilities Xp with 
respect to oscillatory background fields for > 0 and a 
homogeneous background for = 0. The proof of the 
equivalence between Xp and n(p^) is given in appendix A. 
The oscillatory susceptibilities are obtained by evaluat¬ 
ing the appropriate expectation values using noisy esti¬ 
mators, as described in appendix B. Unlike the conven¬ 
tionally used approach, based on position space current- 
current correlators, which mixes information about all 
possible lattice momenta, the present method enables us 
to determine the vacuum polarization with increased pre¬ 
cision for individual low momenta. The low momentum 
region is of relevance for an accurate determination of 
the leading hadronic contribution to the muon anoma¬ 
lous magnetic moment. In principle, the lattice determi¬ 
nation of n(p^) — n(0) at a selected set of low momenta 
can also be combined with experimental results for the 
i?-ratio to increase the accuracy of 

The proposed method not only reduces statistical er¬ 
rors at low momenta but also allows for an indepen¬ 
dent measurement of n(0), instead of having to rely 
on extrapolations of n(p^) from > 0. We discussed 
three different methods to determine the homogeneous 
susceptibility xo = n(0). The most straightforward 
method, which relies only on simulations at zero mag¬ 
netic field (the so-called half-half method), was found 
to suffer from large finite-volume effects of up to 10% 
of the full value. Instead, we combined existing re¬ 
sults on Xo from refs. [48, 58] that are based on sim¬ 
ulations at non-zero background fields. We also tested 
stochastic wall sources to obtain n(0) as the second mo¬ 
ment of a momentum projected current-current correla¬ 
tion function and found that it can compete with the 
accuracy of the homogeneous susceptibility for a suffi¬ 
ciently large number of random sources. It is interesting 
to note that xo can also be obtained via stochastic wall 
sources at finite temperatures, giving direct access to the 
renormalized magnetic susceptibility Xo{T) — Xo(T — 0) 
that enters the QCD equation of state at finite magnetic 
fields [48, 55, 56, 58, 71, 72]. 


10-2 r 


> 

(U 

o 

to 10-^ 
o 
d 
a 




10-5 =■ 


[24] [30] [32] [34] 
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o 


o 

• ; 6 


G 


lattice determinations 


O 


[17] [27] [31] [33] [40] present R- 

work ratio 


FIG. 8. The statistical error of the vacuum polarization at 
low momenta around = 0.03 GeV^ for several lattice stud¬ 
ies in the literature and for the present work (shaded area). 
Open points denote the error of the unsubtracted n(p^), while 
full symbols indicate that of the renormalized nR(p^). Stud¬ 
ies involving only the connected contribution are indicated in 
yellow, while those also taking into account the disconnected 
terms in blue. The determination using the experimental R- 
ratio is also included for comparison (solid green point). 


The method was tested on staggered Nf = 2 -|- 1 flavor 
ensembles with various lattice spacings. Already on a few 
hundred configurations, a statistical accuracy below one 
percent is achieved for n(p^). The disconnected contri¬ 
butions have been included in all cases. Figure 8 shows 
an order-of-magnitude comparison of our statistical ac¬ 
curacy to that of existing calculations in the literature, 
wherever data or figures with error bars are available for 
n at p2 ~ 0.03GeV^ [17, 24, 27, 30-35, 41]. (Note that 
the approach followed in ref. [36] involves parameterizing 
the lattice data for the zero-momentum projected two- 
point function G{t) of equation (9), making a compar¬ 
ison for n difficult.) We remark that this incomplete 
comparison does not distinguish between different lat¬ 
tice volumes, spacings or pion masses but just serves as 
a qualitative indicator of the accuracy. It reveals that 
our statistical errors, obtained on a comparably small 
number of gauge configurations, are by far the smallest 
within the lattice studies shown in figure 8. However, the 
approach of employing the experimental i?-ratio is still 
by about an order of magnitude more accurate. Never¬ 
theless, by applying the methods used in this paper to 
ensembles with substantially higher statistics, the desired 
accuracy may be reached in the near future. 
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Appendix A: Proof of equation (15) 


Below we prove the main result of the paper, equa¬ 
tion (15). We define the free energy density / = 
— log2^/V4, in terms of the partition function Z of the 
system in a four-dimensional volume V 4 . Z is obtained 
evaluating the Euclidean functional integral over the 
gluon, quark and antiquark fields V'/ and ipf, 

r r 

z= , S= (Al) 

J J 


of the action reads 


Sb{B°) = Bj2{0), 


Sb{B^°^’P) = B 
SBiB^'^’^-P) = B 


j 2 {p)-j 2 (-p) 
J2(p) + j2(-p) 


/(2pi), 
/(2ipi), 


(A6) 


where the prime denotes differentiation with respect to 
Pi- Inserting these expressions into the partition func¬ 
tion (Al) and differentiating twice with respect to eB in 
order to obtain the susceptibilities (14) results in 


Xo = ^ (j2i0)j2{0) 


cos ^ ^ 

F 4 dpi 

sin ^ 

V 4 dpi 


hip) - hi-p) 
j 2 ip)+hi-p) 


(A7) 


where the action S is the integral of the Lagrange density 
.if. Without loss of generality, the magnetic field of equa¬ 
tion (13) is chosen to point in the third spatial direction 
and is generated by a vector potential curl A^ = B^: 

A^ = JdxiBP, AP = AP = AP = 0. (A2) 


Note that terms containing the squares of expectation 
values — e.g., (ji (0))^ for xo — vanish due to parity 
symmetry B —B and thus do not appear in equa¬ 
tion (A7). 

Comparing equations (A4) and (A7) shows that 

xr + xf = n(p2). (A8) 


Here the superscript p indicates the oscillation frequency 
of the magnetic field, cf. equation (13). The vector po¬ 
tential enters the Lagrange density via minimal coupling: 


Nf 

f — l 

0/ = 7m + ’i'lfAP^) , 


where .ifg is the gluonic Lagrangian and m/ denote the 
quark masses. 

In equation (A2) we chose a gauge, in which the photon 
vector potential only couples to the second component 
j 2 of the electromagnetic current. Therefore, this back¬ 
ground probes the n 22 (p) entry of the vacuum polariza¬ 
tion tensor, where we orient the momentum p to point in 
the ^-direction: p = (pi,0,0,0). In this case, employing 
equation (1), the vacuum polarization (II) simplifies to 


n22(p) = ^ (^hip)hi-p)) = -pi n(p^) • (A4) 


For reasons that will become clear in a moment, we con¬ 
sider two different oscillatory background fields 


H™’P(x) = B sinipx ), B^°^’P{x) = B cos{px ), 

(A5) 

and denote the corresponding susceptibilities accordingly 
as Xp" and Xp°"- 

Integrating the Lagrange density (A3) and going to 
momentum space, the magnetic field-dependent part Sb 


In the zero momentum limit, the oscillatory magnetic 
fields satisfy 

lim , lim B^'^'^’P = 0, (A9) 

p —>-0 p —>-0 

which, together with equation (A8), implies for the ho¬ 
mogeneous case 


xo = n(o). (Aio) 

Furthermore, the cos- and sin-type magnetic fields only 
differ in a phase and are equivalent due to translational 
invariance. Therefore, the two oscillatory susceptibilities 
coincide, giving: 

2Xp = n(p"). (All) 

Note that the equivalence of x™ and x™^ only holds for 
non-zero momenta and breaks down at p = 0. In addi¬ 
tion, on the periodic lattice the two oscillatory suscep¬ 
tibilities differ at the maximal momentum Pmax = 'xja 
where the cos-type vector potential becomes zero on all 
lattice sites and thus Xp°l„ vanishes identically. (Still, 
equation (A8) holds even at this momentum.) 

Relations (AIO) and (All) represent the basis of our 
analysis to obtain the vacuum polarization function from 
magnetic susceptibilities. We remark that implementing 
equations (A2), (A4) and (A9) may be thought of as us¬ 
ing (5-sources in momentum rather than in position space 
when computing n(p^). A similar idea to relate hadronic 
matrix elements to the response to background fields was 
also discussed in ref. [73]. 
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Appendix B: Implementation of the susceptibilities 


In this appendix we present the details of the lattice 
computation of the susceptibilities (14). First of all we 
have to address the implications of magnetic flux quan¬ 
tization. 


In a finite periodic volume, the magnetic flux <i> 
through the perpendicular plane Li x L 2 (the magnetic 
field is oriented in the third spatial direction) is quan¬ 
tized [74], 


$ = 


J dxi dx 2 eB = dirNf,, 


Nf, G Z , 


(Bl) 


where we exploited that the smallest electric charge in the 
system equals qd = —e/3. Thus, flux quantization pro¬ 
hibits direct differentiation with respect to the amplitude 
of the magnetic field, unless the flux identically vanishes. 
For the oscillatory field BP{x) of equation (13) this is 
indeed the case, making the differentiation with respect 
to B straightforward. For the homogeneous background, 
the flux is non-zero and, thus, B becomes a discrete vari¬ 
able. Various methods to calculate xo summarized 
in refs. [48, 54]. 


After integrating out the quark fields, the lattice par¬ 
tition function becomes an integral over the gluonic links 






r 

z= VUf.e-^oJlidetMPy, + mf . 

(B2) 

Here we employed (rooted) staggered quarks to discretize 
the fermion matrix but the method trivially gener¬ 
alizes to different discretizations. The U(l) vector po¬ 
tential of equations (A2) and (A3) enters via the 
substitution 


U 2 {.xi) ^ U 2 {xx) ■ = U 2 {xi) ■ (B3) 


where in the second step we inserted the vector potential 
for the cos-type magnetic field with momentum pi in 
the first spatial direction. The improvement pi 1 —>■ pi is 
carried out in the denominator of the exponent, similarly 
as in the conventional approach, cf. equation (10). The 
derivative with respect to B is then obtained as 


xp=/-(d+ 


V 4 \ ^ d(eB) / ’ 


where 


f(",') 




(B4) 


(B5) 


and the dot denotes differentiation with respect to the 
combination q/B at B = 0. Having taken the derivative 
at B — 0, we can exploit the equality of the up and down 
quark matrices = Md = due to the coincident 
light quark masses. Then, the susceptibility reads (sup¬ 
pressing the index p and using the electric charge values 
g«/2 = -qd = -Qs = e/3): 


Xp — 




^^tr tr + ^tr tr - ^tr tr , 


(B6) 


where, like in equation (B5), the pre-factors 1/4 and 1/16 
are due to the use of rooted staggered fermions. 

The first expectation value on the right hand side is 
the connected contribution, whereas the second one is 
the disconnected term. The traces are measured via a 
set of noisy estimators 5^, j = 1... A/. Taking into ac¬ 
count the cyclicity of the trace, the total number of nec¬ 
essary inversions is 4A/ (twice for the light and twice for 
the strange quark matrix). For the calculation of A/ dif¬ 
ferent momenta, some of the solutions can be recycled. 
This results in the total number of required inversions 
A^inv = 2Aj(l -I- Np), where the pre-factor 2 again is due 


to M~^ ^ We then have independent estimates 

for the connected contribution. Using different stochas¬ 
tic sources for the strange and for the light quarks, we 
obtain A| estimates of the last disconnected term within 
equation (B6), while for the two non-flavor mixing dis¬ 
connected terms we can only exploit N^{N^ — 1)/2 inde¬ 
pendent variations. 
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Appendix C: Implementation of random wall sources 

Below we specify the details of the calculation of n(p^) 
and of n(0) via stochastic wall sources. In this approach, 
one calculates the current-current correlator in coordi¬ 
nate space, and performs the Fourier transformation sub¬ 
sequently. Care has to be taken in defining the currents 
and especially their product at the same position. Usu¬ 
ally in the literature the conserved current is considered 


and the contact term is subtracted in order for the lattice 
Ward identity (10) to hold [75]. Another possibility is to 
consider the product of conserved and local currents as 
was done in ref. [31]. 

Here we demonstrate how the subtraction of the con¬ 
tact term can be obtained automatically if the current is 
defined using a background U(l) field A^. For simplicity, 
we again consider the fj, = 2 component of the currents 
and take the distance between the insertions to point in 
the first direction. Then we get 


{h{x)j 2 {y)) = 


d'^ log Z 


dA2{x)dA2{y) 


A 2=0 


2 

^1 tr 
e 


M 


-1 


d^M 


/ 


-(y-tr 

16 \ ^ e 


M: 


f d{qfA2{x))d{qfA2{y)) 
dMf 


Sx.if - M 


d{qfA 2 {x))_ 




-1 dMf 
f d{qfA 2 {x)) 

dMf> 


M 


-1 


dMf 


d{qfA 2 {y)) 


(Cl) 


y d{qf,A 2 {y))_ 


Notice that the first term arises due to the fact that the background field enters M/ in the exponential form 
and it only contributes ii x = y. Now we define 


E 


dMf 

d{qfA 2 {x)) 


72 '^xi , 


3 &^Mf 

" d{yfA 2 {x))d{qfA 2 {x)) = 


(C2) 


where Vx^ is the projector on the slice of the lattice where the first spatial coordinate equals xi. Here, 72 is the 
staggered discretization of the second Dirac matrix and T 2 its equivalent with the Hermitian conjugate links multiplied 
by minus one, 


il 2 )xy 

{T 2 )xy 


V2{x)U2{x)6y^^^^^ 


r] 2 ix-a2)U^{x- a2) , 


(C3) 


and r]f 2 denote the staggered phases. With these definitions we obtain for the two-point function (9) — with the 
temporal direction replaced by the first spatial direction: 


G{xi - yi) = ^ 2 E E ih{x)j2{y)) = y E (^) ^j2Vx,Mj: ^-f2Vy, 


X2.X3,X4 2 / 2 , 2 / 4 , 2/4 


(C4) 


-f 


16 


E ^ tr ^72^X1 / ( E ^ Mj:^-f 2 Vy, 


f 


q^ 

e 


All source positions yi can be averaged over, keeping the 
distance xi — yi fixed, to increase statistics. Inserting the 
electric charges and taking into account the degeneracy 
of the light quark masses, this expression can be sim¬ 
plified, in analogy to equation (B6). For its evaluation 
we again use noisy estimators {j = 1 .. .N^) that are 
projected using the V operators. One technical issue is 
the treatment of the second term in the connected contri¬ 
bution of equation (C4). Exploiting the ? 75 -Hermiticity 
Mj = rj^Mfrj^ of the staggered fermion matrix, the fact 
that = V and that the term in question is real, we 


arrive at 

^jT^yil2Mf 'Pxi"t2Mj Vy^^j 

(C5) 

This demonstrates how this term can be obtained for a 
fixed source position yi and any sink position xi using 
only two inversions. One of these inversions can also 
be reused for the calculation of the contact term involv¬ 
ing T 2 and for the traces in the disconnected term of 
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equation (C4). The number of necessary inversions is 
7Vi„v = 4N^. ^ 

Putting all this together, the vacuum polarizations at 
finite and at zero momentum equal (cf. equations (1), (8) 
and (10)), 

X, 

n (C6) 

n(0) = 

ail 

where pi = (2/a) sin(api/2) is the lattice momentum and 


is a quadratic function consistent with the boundary con¬ 
ditions for a periodic lattice with linear size Li. We men¬ 
tion that the separation x — y of equation (C4) is usually 
chosen to lie in the temporal direction, as indicated in 
equation (8). In our setup it points in the first spatial 
direction to make the connection to the magnetic sus¬ 
ceptibilities - involving xi-dependent phases, cf. equa¬ 
tion (B3) - more transparent. 


f{xi) = 



< Ti/2, 

xi)^, otherwise 


(C7) 
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